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1 Introduction 

Kinetic theory can be viewed as an intermediate step between a detailed 
microscopic analysis of a many-body system and the corresponding phe- 
nomenological macroscopic description. In kinetic theory the main objective 
is to derive and solve the kinetic equation for the one-particle distribution 
function, thus obtaining information about the system properties. In the 
context of dilute gases, the Boltzmann equation (BE) [jjj provides the ade- 
quate framework for studying states arbitrarily far from equilibrium. Exact 
solutions to this equation are rare, but a great deal of information can be 
obtained from simplified kinetic models || or from simulation Monte Carlo 
methods ||. On the other hand, the assumptions implicit in the BE are 
only physically justified in the low-density limit. As the density increases, 
structural effects become important, potential contributions to the fluxes 
dominate, and the BE is no longer adequate. There is no general kinetic 
equation valid for finite densities. A singular exception, however, is the 
idealized system of hard spheres of diameter a, for which Enskog proposed 
a semi-phenomenological equation [Q by introducing two crucial changes 
in the Boltzmann collision integral: (a) the centers of two colliding parti- 
cles are separated by a distance equal to c; (b) the collision frequency is 
increased by a factor that accounts for the spatial correlation between the 
two colliding molecules. Although the Enskog equation (EE) also ignores 
the correlations in the velocities before collision (stosszahlansatz) , it leads 
to transport coefficients that are in good agreement with experimental and 
simulation values for a wide range of densities. In addition, the revised 
Enskog theory (RET) Q is asymptotically exact at short times and there- 
fore has no limitations on density or space scale in that limit. Moreover, it 
admits both fluid and crystal equilibrium states as stationary solutions. 
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The mathematical complexity of the EE has hindered practical applications. 
As in the case of BE, two different approaches have been proposed to cope 
with this problem. First, a Monte Carlo algorithm has been introduced to 
solve numerically the EE in the same spirit as the DSMC method of 
solving the BE ||; second, a simple kinetic model that retains the main 
features of the EE has been constructed |6|, ^ ; it reduces in the low density 
limit to the simplest kinetic model of the BE, the Bhatnagar-Gross-Krook 
(BGK) model. Both approaches have demonstrated to succeed in capturing 
the essential properties of the EE and have great potential for a new under- 
standing of nonequilibrium systems under conditions accessible previously 
only by molecular dynamics simulation. 

Two-dimensional systems often serve as prototypes to investigate some 
physical properties also present in real systems. In particular, many of 
the peculiarities of a hard-sphere fluid in far from equilibrium states are 
expected to be present in a hard-disk system with a similar value of the 
packing fraction. Obviously, the calculations usually become easier from 
both theoretical and computational points of view. In this paper we calcu- 
late the rheological properties of a dense fluid of hard disks under shear far 
from equilibrium by using a kinetic model of the EE. The results are com- 
pared with Monte Carlo simulations. As happens in the three-dimensional 
case jjj, the comparison shows an excellent agreement at all densities and 
shear rates considered. 



2 The Enskog Equation for a Hard-Disk Fluid 
under Uniform Shear Flow 

The uniform shear flow (USF) is one of the few inhomogeneous states for 
which exact results can be obtained far from equilibrium, and therefore is of 
great significance in providing insight for the type of phenomena that occur 
under extreme state conditions. Moreover, it has been studied extensively 
by molecular dynamics simulation to analyze rheological properties in simple 
atomic fluids. The macroscopic state is characterized by a constant density 
n, a uniform temperature T, and a linear flow field: u(r) = a-r = ayx, where 
a = axy. The (constant) shear rate a is a single parameter that can be cho- 
sen to drive the system arbitrarily far from equilibrium. The shear produces 
viscous heating that is compensated by an external nonconservative force 
(thermostat), F = — ma(o)V, where m is the mass of a particle, V = v u 
is the peculiar velocity, and the thermostat parameter a is adjusted to as- 
sure that the temperature remains constant. At a microscopic level, the 
USF is characterized by a distribution function, /(r, v, t) — > f(V,t), that 
becomes uniform in the Lagrangian frame of reference. Under the above 
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conditions, the EE for /(V, t) becomes 



where Js[/] is the Enskog collision operator, 

J E [f} = <rx(n) J dVi / d*Q(a-g) (?-g)[/CV / ,t)/(V' 1 ,t)-/(V J *)/CV 1 ,t)] 

(2) 

In the above expression a is the disk diameter, x( n ) is the pair correlation 
function at contact for an equilibrium system with (uniform) density n, <x 
is a unit vector, Q(x) is the Heaviside function, g = V — Vi — era • er, 
V' = V - (<r • g)<x, and = Vi + 2aa ■ er + (<r • g)<x. 
The most relevant transport quantity is the steady-state pressure tensor P, 
which measures shear and normal stresses. It has a kinetic part, P k , and a 
collisional transfer part, P c , that are functionals of / given by 

P k =mj rfVVV/(V) , (3) 
pc = ™fx j dY j j da aa e ( ^.g )(? . g) 2 /(v+aa . ?)/(Vl) . (4) 

Following the standard Chapman-Enskog method [[[] , the Navier-Stokes con- 
stitutive equations are derived and the Newtonian shear viscosity can be 
identified as § 

% , 7T„*.V , 1 _*2_, ^1/2 



ms{n) = — [1 + -^ n *Xj + -^n* x(irmk B T) 1 , (5) 

where 770 = \.Q22(mkBT /n) 1 / 2 /2a is the Boltzmann viscosity, fee is the 
Boltzmann constant, and n* = no 1 is the reduced density. This Navier- 
Stokes shear viscosity is the zero shear rate limit (a — » 0) of a generalized 
transport coefficient r)(n,a) = —P xy /a. Other non-Newtonian effects are 
associated with the differences P xx {n, a)—po(n) and P yy (n, a)—po(n), where 
Po(n) — nksT(\ + \n*x) is the equilibrium hydrostatic pressure. 



3 The Kinetic Model 

Very recently, a kinetic model has been derived by replacing the Enskog 
collision operator with a simpler form that, otherwise, retains the main 
qualitative features. The model has the same domain of applicability and 
preserves the same basic properties (such as local conservation laws and the 
exact equilibrium stationary state for both fluid and crystal phases) as the 
RET. For a detailed account of the kinetic model, we refer the reader to 
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Refs. H |7). In the particular case of the USF, the kinetic model leads to 
the replacement 



JE{f)^-»(f-fl)-fi 



P x V ( m 2 \ 2m 

—a — 1 V X V V A X , 

nk B T \2k B T J k B T v 1 



(6) 



where we have already considered the two-dimensional case. In Eq. v 
represents an effective collision frequency depending on the local density 
and temperature. Here, this parameter is chosen to assure that the low 
density shear viscosity is the same as that from the BE, v = nksTx/vo- I n 
addition, fi(V) is the local equilibrium distribution, P£ y is the xy-element 
of the collisional transfer pressure tensor, as obtained from Eq. ([|), and A xy 
is the collisional moment 

A xy = / dVV x V y J E [f e ] = ~(k B T/m)ina X a (l + ^ , (7) 

where a = ^aa[m / k B T) 1 / 2 . Equation (|l|) together with the substitution 
(|(3]) constitutes now the kinetic equation for the problem. Since the term 
P xy is a functional of /, Eq. ([[]) [along with (||)] is still a highly nonlinear 
integro-differential equation. 

In order to ease the notation, we choose units such that v = 1, m = 1, and 
2k B T/m — 1, and define the dimensionless pressure tensor P* = P/nk B T. 
In these units, a = (v27r/1.022)n*X. Conservation of energy gives the 
thermostat parameter a(a) in terms of P* y (a): 

«(«) = -~P: v (a) ■ (8) 

Taking moments in both sides of the (stationary) kinetic equation, one can 
easily get the kinetic part of the pressure tensor: 

pfc* - 1 I fl2 ~ 2aA *y P k *- l±3SL (2A -a) (9) 

+ (1 + 2a) 2 + a 2 ' *y (1 + 2a) 2 + a 2 1 xy ' ' { ' 

In order to close the mathematical problem, we would need to express P x * 
in terms of a. In principle, this implies to perform the velocity integrals 
(|4[) with the formal solution to the kinetic equation. Instead, we obtain a 
reasonable estimate for P^,* by using the first Sonine approximation 

/ - ft [1 + (v x 2 - v y 2 )(p^ - 1) + 2v x v y p k x ;\ . (10) 

With this approximation, the evaluation of P c * is similar to that of A xy . In 
particular, 

xy 



Pxy = J d v°xV y {{^ + 2a 2 dldl)cYf(aa x d y )-2d x d y P Q 

-(2dl-i?{p k x :-if-Adldlp k x ; 2 ]} 



(11) 
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From Eqs. (||), ([)]) and ( p"l| ) one gets a closed equation for a, whose numerical 
solution can be easily obtained for arbitrary shear rates and densities. 



4 The ESMC Method 

As discussed in the Introduction, a recent method has been developed for 
Monte Carlo simulation of the solution to the EE || . Previous results have 
demonstrated the utility of this Enskog simulation Monte Carlo (ESMC) 
method for studying far from equilibrium states in the regime of low and 
moderate densities. In these conditions, the ESMC method can be even 
more efficient, from a computational point of view, than hard-sphere molec- 
ular dynamics. In addition, the ESMC algorithm reduces to the well-known 
DSMC method j| in the low density limit. 

As applied to the USF, the method proceeds as follows. The one-particle 
distribution function /(V) is represented by the peculiar velocities {V^} of 
a sample of N "simulated" particles. These velocities are updated at integer 
times t = At, 2At, 3At, . . ., where the time-step At is much smaller than 
the mean free time and the inverse shear rate. This is done in two stages: 
free streaming and collisions. The free streaming stage consists of making 

V, e- aAt (Vi - a • ViAt) , (12) 



where the thermostat parameter a is adjusted to assure that the tempera- 
ture, which is computed as T = {m/Nks) , remains constant. In the 
collision stage, a sample of \Nw pairs are chosen at random with equiproba- 
bility, where w is an upper bound estimate of the probability that a particle 
collides in the time interval between t and t + At. For each pair ij belong- 
ing to this sample, the following steps are taken: (1) a given direction <Xy 
is chosen at random with equiprobability; (2) the collision between parti- 
cles i and j is accepted with a probability equal to the ratio Wij/w, where 
Wij = 2iraxnAtQ(a i: j ■ Sij)(a%j ■ gy) and gy = V; - Vj - era ■ ct^-; and 
(3) if the collision is accepted, post-collision velocities are assigned to both 
particles: 

V, -V, • . V , • V , • cT.,<a,, ■ }••,, :• . (13) 

In the case that in one of the collisions Wij > w, the estimate of w is 
updated as w — Wij. In the course of the simulations, the kinetic and 
collisional transfer contributions to the pressure tensor are evaluated. They 
are given as the computational analogs of Eqs. (^) and ([|), i.e. 

N 

P* = 7r£ V « V «. ( 14 ) 

i=l 
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Figure 1: Plot of the normalized shear viscosity Tl/iVo/x) an d its kinetic 
part r] k /(r]o/x) for a = (• • •), a = 0.7 ( — ), and a = 1 ( — ). Lines are from 
the kinetic model; symbols are from Monte Carlo simulations of the EE. 

ran a v-^/— , 

ij 

where the dagger means that the summation is restricted to the accepted 
collisions. Once the steady state is reached, the above quantities are aver- 
aged over time to improve the statistics. 



5 Results and Discussion 

Our objective has been to obtain the pressure tensor P as a function of the 
density and the shear rate by means of the kinetic model, as well as by 
performing Monte Carlo simulations of the EE. 

Figure |l| shows a comparison of the (normalized) non-Newtonian shear vis- 
cosity r\ = —P xy /a as a function of the density parameter n*x for three 
different values of the shear rate. The corresponding kinetic part is also 
shown in the right side of the figure. In all the figures presented in this 
paper, the error bars are smaller than the sizes of the symbols and are not 
drawn. The good agreement indicates that both the kinetic and collisional 
transfer contributions are accurately given by the model. The dotted lines 
correspond to the Navier-Stokes shear viscosity, Eq. (|5|), so that the solid 
and the dashed lines represent non-Newtonian effects. It is quite apparent 
that these effects are much more important at finite density than at zero 
density. It is worthwhile noticing that the shear viscosity presents, in gen- 
eral, a non monotonic behavior as the shear rate increases. More precisely, 
a transition from shear thinning, 77(71,0) < f7Ns( n ) ; to shear thickening, 
77(71,0) > f]Ns(n), takes place when the shear rate is larger than a critical 
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Figure 2: Plot of P* x ( — ), P* y (- -) and pg (• • •) (left side), and the kinetic 
part P x * (right side), for a = 1 . Lines are from the kinetic model; symbols 
are from Monte Carlo simulations of the EE. 



value a c (n) that decreases as the density increases. In particular, a c = 1.0 
at n*x — 2.2, while a c — 0.7 at n*x — 2.7. This transition has been recently 
predicted for a hard-sphere fluid from an analysis of the kinetic model in 
the limit of small shear rates and confirmed by the ESMC method |7| . This 
is an interesting effect not observed for the kinetic part. 
Figure |^ shows the (dimensionless) normal stresses P* x and P* y as functions 
of the density for the shear rate a = 1. The dotted line represents the 
equilibrium hydrostatic pressure p^ = 1 + \n*x of a hard-disk fluid. Again, 
the numerical solution of the model shows an excellent agreement for the 
wide range of densities considered. At high densities the collisional part of 
the pressure tensor dominates and both diagonal elements P xx and P* y tend 
to coincide. In addition, viscometric effects become evident by observing 
the increase of the nonequilibrium "hydrostatic pressure" p* — jtrP* with 
respect to its equilibrium value Pq, a phenomenon usually referred to as 
"shear dilatancy". The kinetic contribution P xx is also plotted in Fig. 2. 
Note that the yy-element can be derived from the consistency condition 
P x * + Py* — 2. We observe that again a good accuracy of the model 
predictions holds for these quantities. 

In this work we have presented results obtained from both a kinetic model 
H and a simulation Monte Carlo method of the EE || for a hard-disk sys- 
tem under shear, arbitrarily far from equilibrium. Although the EE is not 
expected to be accurate for very large densities, we have also considered 
those densities in order to check the kinetic model predictions. Comparison 
between the results obtained from the ESMC method and from the numer- 
ical solution of the model shows an excellent agreement, so that we can 
conclude that the kinetic model gives a good description of the rheological 
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properties throughout the whole shear rate-density plane. In particular, we 
have analyzed non-Newtonian effects beyond Navier-Stokes order by com- 
puting the elements of the pressure tensor. A transition from shear thinning 
to shear thickening is clearly shown by using both approaches. 

Partial support from the DGICYT (Spain) through Grant No. PB97-1501 
and from the Junta de Extremadura (Fondo Social Europeo) through Grant 
No. PRI97C1041 is acknowledged. 
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